      SUBROUTINE GRADEQ (N,MA,A,MB,B,LOW,IGH,CPERM,WK)
CSE
C
CPS        PURPOSE:
C
C             GRADES THE SUBMATRICES OF A AND B GIVEN BY STARTING INDEX
C             LOW AND ENDING INDEX IGH IN THE GENERALIZED EIGENVALUE
C             PROBLEM A*X = (LAMBDA)*B*X BY PERMUTING ROWS AND COLUMNS
C             SUCH THAT THE NORM OF THE I-TH ROW (COLUMN) OF THE A
C             SUBMATRIX DIVIDED BY THE NORM OF THE I-TH ROW (COLUMN)
C             OF THE B SUBMATRIX BECOMES SMALLER AS I INCREASES.
C
CPE
CAS           ***  ARGUMENT LIST  ***
C
C  INPUT PARAMETERS:
C
C  N - THE ORDER OF THE MATRICES A AND B;
C
C  MA - THE ROW (FIRST) DIMENSION OF THE ARRAY A AS SPECIFIED
C       IN THE CALLING PROGRAM;
C
C  A - A REAL TWO-DIMENSIONAL ARRAY WITH ROW DIMENSION MA AND COLUMN
C      DIMENSION AT LEAST N CONTAINING THE MATRIX A AS THE PROBLEM
C      IS DEFINED IN THE FUNCTION STATEMENT ABOVE;
C
C  MB - THE ROW (FIRST) DIMENSION OF THE ARRAY B AS SPECIFIED IN
C       THE CALLING PROGRAM;
C
C  B - A REAL TWO-DIMENSIONAL ARRAY WITH ROW DIMENSION MB AND COLUMN
C      DIMENSION AT LEAST N CONTAINING THE MATRIX B AS THE PROBLEM
C      IS DEFINED IN THE FUNCTION STATEMENT ABOVE;
C
C  LOW - THE INTEGER SPECIFYING THE BEGINNING INDEX OF THE SUBMATRICES
C        TO BE GRADED;
C
C  IGH - THE INTEGER SPECIFYING THE ENDING INDEX OF THE SUBMATRICES
C        TO BE GRADED;
C
C  WK - AN ARRAY OF ANY DIMENSION THAT MUST BE PROVIDED BY THE
C       USER FOR WORK SPACE. WK MUST CONTAIN AT LEAST 2*N LOCATIONS.
C       ONLY LOCATIONS LOW THROUGH IGH AND N+LOW THROUGH N+IGH
C       ARE REFERENCED.
C
C  OUTPUT PARAMETERS:
C
C  A - CONTAINS THE PERMUTED AND GRADED A MATRIX;
C
C  B - CONTAINS THE PERMUTED AND GRADED B MATRIX;
C
C  CPERM - THE REAL ARRAY OF DIMENSION AT LEAST N CONTAINING IN
C          ITS LOW THROUGH IGH LOCATIONS THE COLUMN PERMUTATIONS
C          APPLIED IN GRADING THE SUBMATRICES. THE OTHER LOCATIONS
C          ARE NOT REFERENCED;
C
C  WK - CONTAINS IN ITS LOW THROUGH IGH LOCATIONS THE ROW
C       PERMUTATIONS APPLIED IN GRADING THE SUBMATRICES.
C
CAE
C  REQUIRED SUBPROGRAMS - NONE
C
C  REQUIRED FUNCTIONS - ABS
C
C  REFERENCE - R.C. WARD, BALANCING THE GENERALIZED EIGENVALUE
C              PROBLEM, TECHNICAL REPORT ORNL/CSD-47, OAK RIDGE,
C              TENNESSEE, 1979.
C
C              [SAME TITLE] SIAM J. SCI. STAT. COMPUT., VOL 2, 
C                           NO 2, JUNE 1981, PP 141-152.
C
C ****
C
      DIMENSION A(MA,N),B(MB,N),CPERM(N),WK(N,2)
      DOUBLE PRECISION A,B,CPERM,WK,SUMA,CMAX,RMAX,TEMP,SUMB
      IF (LOW .EQ. IGH) GO TO 505
      IGHM1 = IGH-1
C ****
C     COMPUTE COLUMN NORMS OF A / THOSE OF B
C ****
      DO 420 J=LOW,IGH
         SUMA = 0.D0
         SUMB = 0.D0
         DO 410 I=LOW,IGH
            SUMA = SUMA + DABS(A(I,J))
            SUMB = SUMB + DABS(B(I,J))
  410    CONTINUE
         IF (SUMB .EQ. 0.D0) GO TO 415
         WK(J,2) = SUMA / SUMB
         GO TO 420
  415    WK(J,2) = 1.D37
  420 CONTINUE
C ****
C     PERMUTE COLUMNS TO ORDER THEM BY DECREASING QUOTIENTS
C ****
      DO 450 J=LOW,IGHM1
         CMAX = WK(J,2)
         JM = J
         JP1 = J+1
         DO 430 K=JP1,IGH
            IF (CMAX .GE. WK(K,2)) GO TO 430
            JM = K
            CMAX = WK(K,2)
  430    CONTINUE
         CPERM(J) = JM
         IF (JM .EQ. J) GO TO 450
         TEMP = WK(J,2)
         WK(J,2) = WK(JM,2)
         WK(JM,2) = TEMP
         DO 440 I=1,IGH
            TEMP = B(I,J)
            B(I,J) = B(I,JM)
            B(I,JM) = TEMP
            TEMP = A(I,J)
            A(I,J) = A(I,JM)
            A(I,JM) = TEMP
  440    CONTINUE
  450 CONTINUE
      CPERM(IGH) = IGH
C ****
C     COMPUTE ROW NORMS OF A / THOSE OF B
C ****
      DO 470 I=LOW,IGH
         SUMA = 0.D0
         SUMB = 0.D0
         DO 460 J=LOW,IGH
            SUMA = SUMA + DABS(A(I,J))
            SUMB = SUMB + DABS(B(I,J))
  460    CONTINUE
         IF (SUMB .EQ. 0.D0) GO TO 465
         WK(I,2) = SUMA / SUMB
         GO TO 470
  465    WK(I,2) = 1.D37
  470 CONTINUE
C ****
C     PERMUTE ROWS TO ORDER THEM BY DECREASING QUOTIENTS
C ****
      DO 500 I=LOW,IGHM1
         RMAX = WK(I,2)
         IM = I
         IP1 = I+1
         DO 480 K=IP1,IGH
            IF (RMAX .GE. WK(K,2)) GO TO 480
            IM = K
            RMAX = WK(K,2)
  480    CONTINUE
         WK(I,1) = IM
         IF (IM .EQ. I) GO TO 500
         TEMP = WK(I,2)
         WK(I,2) = WK(IM,2)
         WK(IM,2) = TEMP
         DO 490 J=LOW,N
            TEMP = B(I,J)
            B(I,J) = B(IM,J)
            B(IM,J) = TEMP
            TEMP = A(I,J)
            A(I,J) = A(IM,J)
            A(IM,J) = TEMP
  490    CONTINUE
  500 CONTINUE
      WK(IGH,1) = IGH
  505 CONTINUE
C ****
C     EXIT GRADEQ
C ****
      RETURN
      END
